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ABSTRACT 

We have developed an implicit, multi-group, time-dependent, spherical neutrino 
transport code based on the Feautrier variables, the tangent-ray method, and 
accelerated A iteration. The code achieves high angular resolution, is good to 
0(f/c), is equivalent to a Boltzmann solver (without gravitational redshifts), and 
solves the transport equation at all optical depths with precision. In this paper, 
we present our formulation of the relevant numerics and microphysics and explore 
protoneutron star atmospheres for snapshot post-bounce models. Our major focus is 
on spectra, neutrino-matter heating rates, Eddington factors, angular distributions, 
and phase-space occupancies. In addition, we investigate the influence on neutrino 
spectra and heating of final-state electron blocking, stimulated absorption, velocity 
terms in the transport equation, neutrino-nucleon scattering asymmetry, and weak 
magnetism and recoil effects. Furthermore, we compare the emergent spectra and 
heating rates obtained using full transport with those obtained using representative 
flux-limited transport formulations to gauge their accuracy and viability. Finally, 
we derive useful formulae for the neutrino source strength due to nucleon-nucleon 
bremsstrahlung and determine bremsstrahlung's influence on the emergent and Vr 
neutrino spectra. These studies are in preparation for new calculations of spherically 
symmetric core-collapse supernovae, protoneutron star winds, and neutrino signals. 
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1. Introduction 

With core-collapse supernova explosions, Nature has contrived an elegant means to create 
compact objects, while at the same time seeding the galaxy with the elements of existence. 
Neutrinos play a key role in the phenomena of collapse and explosion, for not only are they 
produced in abundance at the high temperatures and densities achieved in collapse, but they 
arc weakly enough coupled to matter that they transport heat and leptons on a dynamically 
interesting timescale. It is now thought that neutrino heating of the protoneutron star mantle 
drives the supernova explosion (Colgate and White 1966; Bethe and Wilson 1985), but only after 
a post-bounce delay of lOO's of milliseconds to one second. During this delay, the quasi-static 
accreting core, the protoneutron star bounded by the stalled shock wave, radiates neutrinos of all 
species and the net energy deposition in the semi-transparent "gain region" behind the shock plays 
a pivotal role in "igniting" the explosion. However, the precise deposition rate depends upon the 
details of neutrino transfer at low "optical" depths, putting great demands upon the theoretical 
tools employed to calculate the properties of the neutrino radiation fields. The character of that 
radiation depends upon neutrino-matter opacities, neutrino production source terms, and neutrino 
transport. Over the years, neutrino transport theory and the associated microphysics have reached 
a sophisticated level of refinement (Tubbs and Schramm 1975; Lichtenstadt et al. 1978; Bowers 
and Wilson 1982; Maylc, Wilson, and Schramm 1987; Bludman and Schindcr 1988; Bruenn 1985; 
Janka 1991; Mezzacappa and Bruenn 1993a, b; Messer et al. 1998; Yamada, Janka, and Suzuki 
1999). However, despite these efforts, recent progress in modeling supernovae, and new insights 
gained into the character of multi-dimensional neutrino-driven explosions (Her ant et al. 1994; 
Burrows, Hayes, and Pryxell 1995; Janka and Miiller 1996; Mezzacappa et al. 1998), the supernova 
explosion problem is not solved in detail. We know little about the dependence of the ^^Ni yields 
on progentior mass and composition, the iron-peak nucleosynthesis, the explosion energies, the 
nascent pulsar kicks, and the asymmetries and mixing in the explosion debris. Furthermore, we 
still do not know the duration of the post-bounce delay, nor the ensemble of possible neutrino 
signatures. 

In the past, a variety of approximations to the full neutrino transport equations have been 
employed in complex numerical codes meant to simulate stellar collapse and supernova explosions. 
These compromises have been deemed necessary because of the severe CPU constraints of such 
simulations, particularly when those simulations have been multi-dimensional (Herant et al. 1994; 
Burrows, Hayes, and Fryxell 1995; Janka and Miiller 1996). A variety of gray approaches, flux 
limiters, equilibrium assumptions, and approximations to both neutrino source and redistribution 
terms have been employed, sometimes to good effect. However, given the marginality of the 
explosions thus far obtained, the fact that there is as yet no unanimity among theorists concerning 
important issues of principle (cf. Mezzacappa et al. 1998), and the manifest importance of 
neutrinos in collapse phenomenology, a fresh look at neutrino transport and the relevant neutrino 
physics is in order. It is in that spirit that we have constructed an implicit, time-dependent, 
multi-group, multi-angle, multi-species neutrino transfer code to simulate the neutrino radiation 
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fields in stellar collapse and explosion. This code embodies a different computational method 
from that used in the pioneering papers by Bruenn (1985) and Mezzacappa (Mezzacappa and 
Bruenn 1993a, b), but in its use of Feautrier variables and the tangent -ray method it is quite in 
keeping with traditional photon transport and stellar atmospheres simulations (Mihalas 1980; 
Mihalas and Mihalas 1984). In this paper, we describe the basic algorithm, discuss and derive the 
relevant neutrino microphysics, and present high-resolution (in energy, angle, and radius) results 
for representative post-bounce protoneutron star configurations. Hence, for this first paper in our 
series on neutrino transport and microphysics we focus on precision neutrino "atmospheres." We 
present the energy spectra, Eddington factors, angular distributions, phase space densities, and 
neutrino-matter energy couplings. We also derive or discuss the relevant neutrino physics, some 
of it new. We calculate the background neutrino radiation fields for two snapshot models, one of 
which is from the work of Burrows, Hayes, and Fryxell (1995) representing the wind phase that 
follows explosion (Model W), and one of which (our Model BM) was kindly provided to us by Tony 
Mezzacappa and is from a multi-group, flux-limited diffusion simulation by Bruenn (Messer et al. 
1998), 106 milliseconds after the bounce of the core of a Weaver and Woosley (1995) 15 M© star. 
These are meant to exemplify various protoneutron star structures and phases for the purposes of 
a detailed scrutiny of the neutrino sector. Consistent dynamical calculations will follow later in 
the series. 

Neutrinos are the major signatures of the inner turmoil of the dense core of the massive star 
and they carry away the binding energy of the young neutron star, a full 10% of its mass-energy. 
The detection of collapse neutrinos, their "light curve" and spectra, will allow us to follow in real 
time the phenomena of stellar death and birth. The supernova, SN1987A, provided a glimpse of 
what might be possible, but it yielded only 19 events; we can expect the current generation of 
underground neutrino telescopes to collect thousands of events from a galactic supernova. 

In §^, we present the equations and physics of neutrino transport. In §^ we describe our 
implementation of the Feautrier and tangent-ray schemes, and follow this in §^ with a discussion 
of accelerated A iteration and our approach to the implicit coupling of matter with neutrino 
radiation. Section 5 contains a physical derivation of stimulated absorption and §^ summarizes 
the cross sections and source terms we employ for this study. We provide in §0 a derivation of 
the single and pair neutrino rates and spectra due to nucleon-nucleon bremsstrahlung, a process 
that can compete with pair annihilation as a source for i/^, v^, Ur, and Ur neutrinos and that 
to date has not been incorporated into supernova codes. (The consequences of bremsstrahlung 
for the emergent spectra are presented in ^10|.) In §|8|, we present for generic protoneutron 
star configurations our basic results vis a vis emergent energy spectra, luminosities, and energy 
deposition rates (including that due to ui' annihilation). We also explore the dependence of the 
emergent spectra and neutrino heating rates on blocking factors, weak magnetism and recoil, 
aberration and Doppler terms, and stimulated absorption. These terms/effects are frequently 
dropped in simpler schemes. In §P, we highlight the angular dependence of the radiation fields 
and the conceptual limitations of flux limiters that ignore the angular dimension. Moreover, we 
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compare the emergent spectra and heating rates obtained using our full transport code with those 
obtained using representative flux-limiter closures in order to gauge the errors of such approximate 
schemes. 

This paper contains a description of neutrino transfer, our numerical approach, and the new 
results that flow from it. It is also meant to summarize various useful formulae that others, as 
they approach the study of supernova neutrino radiation fields, might employ. In assembling the 
rates and cross sections, we have borrowed from the investigations of Tubbs and Schramm (1975), 
Bruenn (1985), Janka (1991), Mezzacappa and Bruenn (1993a,b,c), Schinder and Shapiro (1982), 
and Bowers and Wilson (1982), but take full responsibility when we have chosen to deviate from 
the literature. 



2. Neutrino Transport Equations 

We have constructed a radiation/hydrodynamic code that solves the three equations of 
hydrodynamics with the equations of multi-group radiative transfer and composition. The hydro 
code is a one-dimensional Lagrangean realization of the explicit Piecewise-Parabolic Method 
(PPM) of Colella k. Woodward (1984) (Fryxell et al. 1991) that is automatically conservative, 
second-order accurate in space and time, and employs a Riemann solver to handle shocks. 
Radiation is coupled to the matter between the hydro updates in an implicit, operator-split 
fashion, employing accelerated A iteration (ALI) to facilitate the convergence both of the transport 
solution and of the temperature and composition changes due to transport. Since in this paper 
we focus on the transport sector of the code and on precision neutrino atmospheres, we postpone 
to a later paper a discussion of the full hydrodynamic technique and of time-dependent results 
in the stellar collapse and supernova context. Here we describe the radiation equations solved, 
the algorithm developed to solve them, and the philosophy behind our methods. In later sections, 
we explore the nature of the neutrino radiation fields in the post-bounce and protoneutron star 
contexts. In addition, we study the infiuence of various terms and physics on the emergent 
neutrino spectra and on the neutrino-matter coupling in the semi-transparent region between the 
neutrinospheres and the stalled shock. Energy deposition in this region is thought to be important 
in igniting and driving the supernova explosion. 

Neutrino transport is not an esoteric subject apart from traditional radiative transfer. The 
same techniques developed for one particle type can be employed for another. For all particles, the 
solution to the Boltzmann equation is sought. What distinguishes neutrino transport and transfer 
arc the number of neutrino species (six), the particular microphysics of the neutrino-matter 
interaction [i.e., cross sections, sources), the Fermi statistics of the neutrinos (manifest only in 
the collision term), and the fact that there is in principle a conserved lepton number. Neutrino 
oscillations can alter this, but given the particular neutrino masses and oscillation angles suggested 
by the recent solar and atmospheric neutrino data (Suzuki 1998; Totsuka et al. 1998), oscillations 
might not dramatically affect supernova dynamics or the neutrino fields in the core (Fuller 
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et al. 1992). (It should be borne in mind that oscillations in the outer envelope of the progenitor 
massive star or between the supernova and the Earth may alter the signal detected in underground 
neutrino telescopes.) 

There are a variety of ways of writing the transport equation for the specific intensity (I^) 
of the radiation field (Mihalas 1980; Mihalas and Mihalas 1984). In principle, the Boltzmann 
equation and the transport equation are equivalent, though the former is written in terms of 
the invariant phase space density (Tu), related one-to-one to the specific intensity through the 
identity: 



(1) 



where g = 1 for massless neutrinos, g = 2 for photons, e is the particle's energy, and the other 
symbols have their standard meanings. Sometimes it is said that the Boltzmann equation is 
more general than the transport equation because it contains a p term that for massless particles 
corresponds to gravitational redshifts. However, there is no reason to exclude such a term from 
the transport equation and we will not engage in such distinctions. 

One form of the transport equation for Ii, in the comoving frame in spherically-symmetric 
geometry is: 
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a is the matter acceleration, /i = cos 6, e is the neutrino energy, and tJi, is the emissivity of the 
medium. The subscript u indicates neutrino energy dependence and ^ is an angular phase function 
for neutrino scattering into the beam. This equation, good to 0{v/c) (where v is the matter 
velocity and c is the speed of light), assumes azimuthal symmetry and contains the appropriate 
redshift, aberration, and advection terms due to matter motion, angular redistribution due to 
scattering into the beam, scattering and absorption out of the beam, and source terms. Eq. (H) 
does not include energy redistribution upon scattering, to be incorporated in a later version of 
the code. The various terms represent the additions and subtractions from the beam, the entire 
equation representing conservation of energy and number. The microphysics and collision terms 
reside on the right-hand-side and the geometry, aberration, advection, and Doppler shift terms 
reside on the left. 

While eq. (|^ contains the relevant terms to 0{v/c), it is a bit awkward to difference. It 
is also a bit ugly and its various terms are not so cleanly distinguished by their physical roles. 
Dropping the acceleration term, we follow Eastman and Pinto (1993) and derive the form of the 
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transport equation we employ in this study. The equation, physicahy equivalent to the Boltzmann 
equation (ignoring gravitational redshifts and the acceleration term) for an individual neutrino 
species, is: 



c Dt dr r d/i r \ ^ \ dine, 

= vu-xuiu + ^J^{^,n')iu{n')dn', (4) 

where Q = dlnv/dlnr — 1 and all other symbols have their standard meanings. $ is a phase 
function for neutrino scattering into the beam. Xu is the total extinction coefficient {= Ka + Hs): 
where Ka and Kg contain contributions from all absorption and scattering processes, respectively: 

Kg = ^ rii (t| and Ka = ^ rii . (5) 

i i 

Eq. (0) can be rewritten as the corresponding Boltzmann equation for T^: 
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which can be mapped directly, term-by-term, into the Boltzmann equation employed by Messer 
et al. (1998) in their recent work on Boltzmann neutrino transfer. Eq. (|6|) is the most useful form 
of the transport equation when studying it using the method of characteristics. 

For neutrinos, the phase function for a scattering process i is well approximated (except for 
iy-e~ scattering) by 

$,(n, n') = $,(n ■n') = {i + 5ift- n') = (i + 5^^,) , (7) 

where 6i is a constant specific to each scattering process and here fi is the angle between the 
incident and outgoing neutrinos. Hence, we can write the differential cross section for a scattering 
process i in terms of the total scattering cross section: 

Subsequently, we drop the superscript s. Eq. (0) implies that the angular redistribution term 
in eq. (Q) becomes 

^^sJu + ■ , (9) 

47r 

where 

5t = . (10) 
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Fy in eq. @ is the neutrino flux and Jy is the zeroth moment defined by 

Ju = \f_^Iudti = -^E,, (11) 

where Ejy is the neutrino energy density. 

Integrating eq. and Q, x eq. over dfl yields the zeroth and first moment equations, 
respectively: 

^ ^^"^ + ^^{r^H,) - —{^Pu -Ju) + -U- (J. + QPu) = <{B. - Ju) (12) 
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K* + Kg - ;^^s'^T^ H„ = - (k* + Ktr) i?^ , (13) 

where H^, P^, and N^, are the first, second, and third angular moments given by 

H, = ^l\l,dfi = ^F, , (14) 
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-Bjy is the equilibrium (black body) spectral energy density times ^ and k* includes the correction 
for stimulated absorption (see Ktr in eq. (13) is the total transport extinction coefficient and 
is defined in terms of the individual transport cross sections as Hfr = X^i nia^ ■ For a particular 
scattering process i, 

= l^{l-li)dn = aAl- . (17) 



dn 

Note that 6p and 5„, the (5jS for neutrino-nucleon scattering, are negative {5p ~ —0.2 and 
6n ~ —0.1) and, hence, that these processes are backward-peaked. The fact that 6p and (5„ are 
negative and, as a consequence, that aj^ is greater than cjj increases the neutrino-matter energy 
coupling rate for a given neutrino flux in the semi-transparent region. This increase in inverse flux 
factor {Jy/Hy) is but one effect that can be studied with the transport tools we have developed 
and are developing. 



Integrating eq. ( |l2|) over energy, we obtain the neutrino energy equation: 
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where E and F are the integrated neutrino energy density and flux, respectively, p is the 
energy-integrated Eddington factor, where pi, = P^/Ju- The sums for all neutrino species 
of the negative of the right-hand-side of eq. (^) and the negative of the e^, integral of the 
right -hand-side of eq. ( p!3[ ) are the energy and momentum source terms in the matter equations. 
The two equations for the rate of change of the electron fraction (Ye) due to e~ / e'^ / Ve / i^e capture 
are: 

pMa-^ =±47r/ <(i?,-J,)±-, (19) 

Vt ± Jo £ 

where the — sign is for the equation and the + sign is for the equation. Eqs. (^), (p^), and 
( Jl^ for each neutrino species, along with eqs. (p^), are the basic neutrino transport equations 
that we solve, p is the mass density and is Avogadro's number. 



3. Method of Solution: Feautrier and Tangent Ray Algorithm 



We solve the moment eqs. (12) and ( fi3D implicitly for Jy and Hy by backwards differencing 
in time the quantities Ju/ p'^^^ and H^/ p^/^ , backwards differencing in hie (according to the slope 
of the characteristic), and combining the spatially diff'erenced equations into one equation for Jy 
which is 2nd-order accurate in r. Standard matrix inversion techniques are employed to obtain 
Jjy, from which Hy is derived using eq. (|l^). This equation is manifestly Lagrangean and by 
solving it the advective derivative is included automatically. Jyj p^l"^ and Jiyj j?!"^ are the natural 
combinations for adiabatic compression or expansion of a relativistic gas. 



Since eqs. (|l^) and (|13|) contain the higher-order angular moments Py and closure 
relations are needed. These are obtained from a formal integration of the full transport equation 
(Q), written in terms of the Feautrier variables: 

[/,(/.) = ^(/,(/.)+/,(-^)) and v;(M) = ^(/.(/i)-/.(-M)) , (20) 

where p ranges from to 1. 

In the isotropic scattering limit [bx = 0), these equations for Uy and Vy are: 

IDUy dVy 1-P\ ^ dUy , ^ ^\<^^n^\TT TT^ J 

TTT + I^Q^^-^ h - 3 - — (1 + Qp^)Uy =riy- XyUy + KgJy (21) 

c Ut or r op r \ ome / 

and 

I DVy dUy 1-P\ ^ dVy , (5 ^ ^ ^ H 2 N T/ T/ ( 00\ 

From the solution of eqs. (H) and (H), we obtain the full radiation field and the higher-order 
moments that are then used in eqs. (|l2|) and ([T3| ) for Jy and Hy. Since eqs. (^) and ( p2|) require 
the lower-order moment Jy (and in principle Hy, cf. eq. P), we iterate this system until we obtain 
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a converged and consistent global solution. Simultaneously, we calculate the A operator that maps 
Sy, the source function, onto J^, and employ accelerated A iteration (Cannon 1973a, b; Scharmer 
1981; Olson, Auer, and Buchler 1986; Eastman and Pinto 1993) to speed the convergence of the 
temperature and composition updates. Independent of the total optical depth, this generally 
requires no more than 2 to 3 iterations to obtain an accuracy of a part in 10^. To maintain 
stability and reflect the density character of Uu{fi) and the flux character of Vy{n)^ we stagger the 
[/[/(/i) and V^(/i) meshes with respect to one another by half a zone. 

It may seem that by solving the moment equations separately and iterating with the solution 



of the transport equation itself and by not focusing simply on the solution of eqs. (21) and 
(|2^ or eq. (^) that we are doing more than is necessary. An advantage of solving the moment 
equations is that they can be differenced to automatically conserve energy. However, enforcing 
energy conservation by construction does not guarantee that the solution obtained is the correct 
one. In fact, it is standard with many "non-conservative" hydrodynamic codes that do not 
difference the equations to conserve energy by construction to employ the degree to which 
energy is in fact conserved with time to assess their accuracy. We have chosen to retain the 
automatic energy-conserving feature, but in dynamical calculations we use electron lepton number 
conservation the global code check. By not differencing the equations to automatically conserve 
lepton number, since the equations we difference certainly do conserve lepton number, this is a 
useful approach. This is akin to employing entropy conservation in adiabatic flow as a check of a 
hydro code that is forced to conserve energy automatically by the particular differencing scheme 
employed. However, for the stationary atmospheres calculations we present here, for which the 
time derviative is zero, this is a moot point. 



The advantage of calculating Ui, and instead of is that equations ( |2l| ) and (^) can 
be differenced in such a way that Vy will go accurately to SfidS^/dTu in the large-optical-depth, 
diffusion limit, and still remain accurate in the optically-thin, free-streaming limit. Schemes based 
directly on eq. (^) or (^ have the correct large-optical-depth behavior for Uu, i.e. 17^ = 3^, but 
have round-off trouble computing V^, which is important if the only estimate of the flux comes 
from integrating /xV^ over angle (Larsen, Morel, and Miller 1987). Typically in such codes above a 
certain optical depth the diffusion limit is assumed and V^, is set equal to ?>^dSy/dTy. 

Numerical methods which solve for the spatial variation of a specific-intensity-like-variable 
{e.g., I[/i, r, t]), such as all discrete ordinate transport methods, suffer from the problem that at 
large optical depth, the flux is not well determined (Morel, Wareing, and Smith 1996). This is a 
well-known problem, and the principle motivation for switching to Feautrier variables and to full 
range moment methods, in which J and H are solved for directly. Any S^r method has a spatial 
truncation error which is proportional to (Ar)'^, where k depends on the spatial discretization 
scheme. As Ar grows, the error in the flux, which is proportional to (Ar)'^, also grows. For 
Feautrier variables, on the other hand, finite difference systems of equations can be derived which 
are at least 2'nd-order accurate in r. These give an accurate representation of the radiation field 
in the free streaming limit and go naturally over to the diffusion limit when At is large. 
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To solve eqs. (^) and (|2^), we employ the tangent-ray method (Schinder and Bludman 
1989; Mihalas and Mihalas 1984). At a reference radial zone, tangents are constructed to each 
of the interior zones. The angles of the tangent rays to the normal at the reference zone define 
the angular grid at that zone on which the angular integrations are performed. Eqs. (^) and 
( p^ are integrated along each tangent ray. If there are nx radial zones, the radiation field at the 
outer zone is resolved with nx — 1 angles; as you move inward the number of angles employed 
decreases linearly. Hence, if there are 100 radial zones, there are as many as 99 angular bins. With 
reasonable radial gridding, this approach can provide exquisite angular resolution, particularly 
for forward-peaked radiation fields, but at the cost of increased computational overhead. For 
instance, we have tested our implementation of the tangent -ray method with the Kosirev (1934) 
(spherical Milne) problem for which the absorptive opacity is assumed to be a power law in radius 
(k = 1/r"). For a variety of integer power laws {e.g., n = 1.1, 1.3, 1.5, 2, 3, 4), with from 100 to 500 
radial zones, the tangent-ray method is superior to many implementations of the discrete ordinate 
method (Schinder and Bludman 1989). However, care must be taken to avoid purely geometrical 
zoning, for which r„+i = ar„, since such zoning biases the angular binning in a systematic way. 
The result can be that the Eddington and fiux factors asymptote at infinity to between 0.96 and 
0.98, and not to 1.0. However, purely geometrical zoning is easily avoided in real calculations. 
Note that the increase in angular resolution with radius which comes naturally from this procedure 
is quite appropriate to spherical symmetry. Note also that the tangent-ray method is in some sense 
automatically adaptive for moving grids. At r = 0, the radiation field is by definition symmetric 
and needs no angular resolution. As r becomes large, a small bright central source is increasingly 
finely resolved. 

For dynamical, time-dependent calculations, solving eqs. (^) and (^) by the tangent-ray 
method at each timestep can be time-consuming, but manageable. Fortunately, as long as the 
2'nd and 3'rd angular moments do not change quickly, one need not solve eqs. (|2l|) and ( |2^ ) at 



every timestep. Frequently, the solution to the moment equations (12) and (|T3|) with previous 
values of pi, and (= N^/H^) can be quite accurate. In fact, this is most often the case, since 
the neutrino radiation fields rarely change on timescales shorter than ~0.1 milliseconds, whereas, 
due to the explicit nature of the hydro portion of the code, the Courant timesteps are often near 
one microsecond. Hence, during much of the pre-explosion delay and protoneutron star phases, 
as well as during much of the core collapse phase, it is quite legitimate to solve for the Eddington 
factors only every few steps. An exception is during very dynamical phases such as shock breakout 
through the neutrinospheres. 

Currently, there are two approximations in our algorithm for solving the transport equations 
( |2l|) and (^). The first is that in calculating the radiation angular moments we assume that 
scattering into the beam is isotropic, while maintaining the correct transport cross section in 
the Hi, moment eq. (p!3|). To approximately incorporate the effects of anisotropic scattering, we 
employ the transport cross sections, as described and discussed above. Straightforward methods 
for solving the fully anisotropic problem using the Feautrier variables and the approaches outlined 
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in this paper will be described elsewhere. The second approximation involves the assumption of 
linearity of the characteristics. The 2'nd term in both eqs. ( |2l| ) and ( |2^ ) which is proportional 
to d/djji comes from the aberration experienced in going from the local to a nearby rest frame. 
The characteristics are not perfectly straight, which can make the calculation more difficult. 
One cannot simply integrate along a straight -line impact ray. However, these terms are often 
insignificant because we require only an estimate of Up and Vp and are using them to compute only 
the closure coefficients, and g^. Rather than just ignore these two terms, we have substituted 

- ^/?gf/. (23) 



5/x 

in equation (|2^) and 



l^eQ,?YL ^ ^-l^^Q.V^ (24) 

r a/x r 

in equation (p2[). The substitution in eq. (|2^) is derived by integrating the left -hand-side by parts 
and the substitution in eq. (^) is derived by integrating ^ times the left-hand-side by parts. 
Importantly, the two terms in ( ]23| ) integrate to the same thing and /i times the two terms in eq. 
( p^ ) integrate to the same thing. Therefore, both modified equations reduce to the energy and 
momentum conservation equations. 

In sum, we solve two coupled moment equations for the mean intensity and flux of the 
radiation field. These are the fundamental results of the transport calculation. They are solved 
by an Eddington factor iteration wherein a set of angle-dependent equations consistent with the 
moment equations are solved for the intensity given a constant source function, and this intensity 
is used to determine the closure factors in the moment solution. 



4. Implicit Coupling to Matter and Accelerated A Iteration 

Though we are not highlighting in this paper time-dependent calculations, we think it 
useful to include a discussion of the technique we employ to couple matter with neutrinos. 
This is done implicitly in operator-split fashion, after each hydro update. For each neutrino 
species, the scattering and absorption opacities and the emissivities are calculated and fed into 
the transport solver. A fully converged solution of the transport equations is obtained and this 
is used to calculate the various terms needed for the implicit update of the temperature and 
Ye, due to transport. In particular, the derivatives with respect to temperature and of the 
right-hand-sides of eqs. (|T8| ) and ([l9t ) are calculated. For the implicit temperature update at each 
radial zone, i, a backward-differenced matter energy equation like: 

pCv'-^ K^^ = -^"/o iv'-<J'^)de-^nAT,l {%--^J^-<if)de (25) 

is constructed, where AT, is the temperature change between iterations, t!^^"^ is the new 
temperature, is the old temperature, Cy is the specific heat, r]' is not corrected for final-state 
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neutrino blocking (§^), and At is the timestep. (In fact, the matter energy is a function of both T 
and Yf. and there is an extra term in the temperature update equation to account for the entropy 
change due to the Y^, composition change. That term has been dropped here for clarity, but not in 
the computations.) 

The subtlety with eq. (p5|) lies in the term. In general, Jj equals ^ijSj at each frequency 
or energy, where the A operator is a matrix coupling different zones. Hence, eq. (25) is a matrix 
equation with 

dJi . dSj , , 



For simplicity in eq. (|25|) , we have dropped in eq. ( P6[) the T operator that couples energy groups. 
Though we have the option in the code of calculating the full A matrix, we use only the diagonal 
and the two adjacent off-diagonals. It is this truncated tridiagonal A operator that we actually 
employ. 

Since the Sy we use in eq. (p6|), perhaps a bit idiosyncratically, equals rj' /k* , is given by 



/<, (27) 



and this is employed to derive: 



dJi dSiJ^i dSi^i dSi 

— =Ai,i+i +Ai,i_i +Ai,j— . (28) 



Plugging eq. (28) into eq. (^), we solve for ATj by inverting the tridiagonal matrix. For the Uf, 
and i>e species, a similar procedure is followed to obtain AY^ from eq. (|l9|). Note that the integral 
over neutrino energy is performed before the T and Yg updates, which are not attempted for each 
energy group individually. 

Once ATj and AY^ are obtained, T^+^ is set equal to Tj' + ATj and y^^^+i'* is set equal to 
ykyi _|_ ^Y^. We then loop back to obtain a new transport solution with the new temperature and 
Ye- This procedure is iterated until ATi/Ti and AY^/Y^ are suitably small (normally 10^^) for all 
zones, at which time we are left with a completely consistent set of Jy, H^, Uu, Vy, T, and Y^ . T-' 
and Y^ are not changed during the iteration. The total number of iterations varies between 1 and 
7, the latter only when Y^ is changing fast in either the or the modules. The various neutrino 
fluids are updated in series, not in parallel and we generally follow three species: fg, Ve, and "f^," 
the latter representing the sum of f^, i>^, Vr, and T>r neutrinos. Bunching these four neutrino 
species into one assumes that their cross sections and source terms are identical, which technically 
is false, but quantitatively reasonable). Note that to achieve stable iteration, it is essential for the 



derivatives in eq. (25) to be accurate. Among other things, this requires good derviatives of fi 
(= l^n ~ f^p) with respect to T and Y^. Analytic derivatives are preferred, but numerical derivatives 
for most quantities seem to work. 
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As stated previously, to achieve rapid convergence of the transport iteration we employ ALT, 
accelerated A iteration (Cannon 1973a,b; Scharmer 1981; Olson, Auer, and Buchler 1986). This 
entails an approximation that improves during the iteration. In particular, we use 

4+^ = 4 + A^{S^+^ - S^) , (29) 

where A'^ is the retarded A operator. We use only its diagonal and first off-diagonal terms. 
This procedure accelerates and stabilizes the iteration, even if the optical depth is large and the 
scattering albedo is high. Note that one cannot iterate on the full A matrix (the inverse of the 
matrix representation of the finite-difference transport equations) because its eigenvalues are very 
close to the unit circle and the iteration stabilizes instead of converging. Subtracting off a piece of 
the A matrix and lagging the iteration of that piece allows the iteration to converge much more 
rapidly. 



5. Stimulated Absorption 



The concept of stimulated emission for photons is well understood and studied, but the 
corresponding concept of stimulated absorption for neutrinos is not so well appreciated. This may 
be because its simple origin in Fermi blocking and the Pauli exclusion principle in the context 
of net emission is not often explained. The net emission of a neutrino is simply the difference 
between the emissivity and the absorption of the medium: 



(30) 



All absorption processes involving fermions will be inhibited by Pauli blocking due to final-state 

occupancy. Hence, rjy in eqs. (|30| ) and (§) includes a blocking term, (1 — !Fu) (Bruenn 1985). 

is the invariant distribution function for the neutrino, whether or not it is in chemical equilibrium. 

We can derive stimulated absorption using Fermi's Golden rule. For example, the net collision 
term for the process, e~p^ is: 



(27r)3 J (27r)3 J (27r) 



d?Pu^ f d^Pn f d^pp f d?p, 



(31) 



X r.[Uen ^ e p) (27r) 6 (p^.^ + Pn - Pp - Pe) , 

where p is a four-vector and 

E{v,n ^ e~p) = Tu,Tn{l - - Tp) - TeTp{l - Tn){l - T^:) . (32) 

The final-state blocking terms in eq. ( [3^ ) are manifest, in particular that for the v^. neutrino. 
Algebraic manipulations convert H(fen <-> e~p) in eq. (32) into: 

T' 



E{uen ^ e p) = Tn{l - J^e){^ - J^p) 
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{K^-^v^ , (33) 



where 

= [e(^^e-(/^.-A))/3 ^ -^pi ^34^ 

is an equilibrium distribution function for the neutrino and it has been assumed that only the 
electron, proton, and neutron are in thermal equilibrium. Note that in T'y^ there is no explicit 
reference to a neutrino chemical potential, though of course in beta equilibrium it is equal to 
\ie — A- There is no need to construct or refer to a neutrino chemical potential in neutrino transfer. 

Using eq. (|l]), we see that eq. (|^) naturally leads to: 

Jnet = - I,) = kI{B, - . (35) 



This is akin to the right-hand-side of eq. (12). If neutrinos were bosons, we would have found 
a {\ + J-l) in the denominator, but the form of eq. (^) in which Jjy is manifestly driven to 
By, the equilibrium intensity, would have been retained. From eqs. ( |33|) and (|35|), we see that 
the stimulated absorption correction to is 1/(1 — J-^). If we want to retain the form of the 
collision term as expressed in eqs. ( pO|) or (Q), the physics is unaltered and stimulated absorption 
is not needed as a concept, as long as rjy in eq. (^) contains the neutrino Pauli blocking term, 
(1 — J-y), without the prime. However, by writing the collision term in the form of eq. (35), 



with Ka corrected for stimulated absorption, we have a net source term that clearly drives 
to equilibrium. The timescale is 1/ck* . Though the derivation of the stimulated absorption 
correction we have provided here is for the <-> e^p process, this correction is quite general and 
applies to all neutrino absorption opacities. 

Kirchhoff 's Law, expressing detailed balance, is: 

K'a = Vu/Bu or K* = 'q'y/By , (36) 

where r/^ is not corrected for final-state neutrino blocking. Furthermore, the net emissivity can be 
written as the sum of its spontaneous and induced components: 



where -|- or — is used for bosons or fermions, respectively. 



(37) 



6. Neutrino Cross Sections 

Neutrino-matter cross sections, both for scattering and for absorption, play the central role 
in neutrino transport. The major processes are the super-allowed charged-current absorptions 
of Ue and z^g neutrinos on free nucleons, neutral-current scattering off of free nucleons (Schinder 
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1990; Janka et al. 1996; Burrows and Sawyer 1998ab; Reddy et al. 1998; Yamada 1998), alpha 
particles, and nuclei (Freedman 1974; Leinson et al. 1988; Horowitz 1997; Burrows et al. 1981), 
neutrino-electron/positron scattering (Schinder and Shapiro 1982ab,1983), neutrino-nucleus 
absorption, neutrino-neutrino scattering, neutrino-antineutrino absorption (Janka 1991), and 
the inverses of various neutrino production processes such as nucleon-nucleon bremsstrahlung 
and the modified URCA process {v^ + n + n — > e~ + p + n). Compared with photon-matter 
interactions, neutrino-matter interactions are relatively simple functions of incident neutrino 
energy. Resonances play little or no role and continuum processes dominate. Nice summaries of the 
various neutrino cross sections of relevance in supernova theory are given in Tubbs and Schramm 
(1975) and in Bruenn (1985). In particular, Bruenn (1985) discusses in detail neutrino-electron 
scattering and neutrino-antineutrino processes (see also Dicus 1972) using the full energy 
redistribution formalism. He also provides a serviceable approximation to the neutrino-nucleus 
absorption cross section (Fuller 1982; Fuller, Fowler, and Newman 1983; Aufderheide et al. 1994). 
Recall that for a neutrino energy of ~10 MeV the ratio of the charged-current cross section to the 
fe -electron scattering cross section is ~100. However, neutrino-electron scattering does play a 
role, along with neutrino-nucleon scattering and nucleon-nucleon bremsstrahlung, in the energy 
equilibration of emergent neutrinos, though the relative contribution of each has yet to be 
determined. In this context, our current lack of an energy redistribution algorithm should be borne 
in mind. Nevertheless, our general conclusions in §10 concerning the neutrinos, their softer 
than previously-believed spectra, the likely role of bremsstrahlung in their production, and the 
consequences of their high scattering albedos, will only be strengthened when competent energy 
redistribution is included. 



6.1. Charged-Current Absorption 



The cross section per baryon for either fg or Vf, absorption on free nucleons is larger than that 
for any other process. Given the large abundances of free neutrons and protons in protoneutron 
star atmospheres, these processes are central to Vf, neutrino transport. A convenient reference 
neutrino cross section is cTq! given by 

4G2(mec2^2 



0"n 



~ 1.705 X 10"^"^ cm^ 



The total fe + n — > e~ + p absorption cross section is then given by 



0"o 



1 + 35^ 



e^. + A 



np 



~\~ ^np 



1/2 



Wm 



(38) 



(39) 



The corresponding absorption cross section for the process, + P 



e"*" + n, is 



l + 3ffi\ (e 



A 



np 



A 



np 



1/2 



w, 



M 



(40) 
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QA is the axial-vector coupling constant (~ —1.26) and Anp = rrinC^ — nipC^ = 1.29332 MeV. Wm 
is the correction for weak magnetism and recoil (Vogel 1984), never before included in supernova 
simulations, and is approximately equal to (1 + l.le,^^ /ninC^) for i/g absorption on neutrons. At 
£i,^ = 20 MeV, this correction is only ~ 2.5%. The corresponding correction (Wj^j) for i>e neutrino 
absorption on protons is (1 — 7.1epg/m„c^), which at 20 MeV is a large —15%. To calculate k*, 
o"°^„ and (Tp^p must be multiplied by the appropriate stimulated absorption correction, 1/(1 — J^IJ 
or 1/(1 — J^'p^)- Furthermore, final-state blocking by either electrons or positrons and either 
protons or neutrons (a la eq. [3^ ) must be included. The consequences of these various terms 
for the neutrino spectra and neutrino-matter heating rates are explored in §|8|. Note that the 
sign of yUe — /i in the stimulated absorption correction for i^g neutrinos is flipped, as is the sign 
of He in the positron blocking term. Note also that the + P — ^ e"*" + n process dominates 
the supernova neutrino signal in proton-rich underground neutrino telescopes on Earth, such as 
Super Kamiokande, LVD, and MACRO, a fact that emphasizes the interesting complementarities 
between emission at the supernova and detection in Cerenkov and scintillation facilities. 

7. Nucleon— Nucleon Bremsstrahlung 

A production process for neutrino/anti-neutrino pairs that has received but little 
attention to date in the supernova context is neutral-current nucleon-nucleon bremsstrahlung 
(ni + n2 ^ + n4 uU). Its importance in the cooling of old neutron stars, for which the 
nucleons are quite degenerate, has been recognized for years (Flowers, Sutherland, and Bond 
1975), but only in the last few years has it been studied for its potential importance in the 
atmospheres of protoneutron stars and supernovae (Suzuki 1993; Burrows 1997; Hannestad 
and Raffelt 1998). As a consequence, it has never before been incorporated into supernova 
codes. Neutron-neutron, proton-proton, and neutron-proton bremsstrahlung are all important, 
with the latter the most important for symmetric matter. As a source of and i>e neutrinos, 
nucleon-nucleon bremsstrahlung can not compete with the charged-current capture processes. 
However, for a range of temperatures and densities realized in supernova cores, it may compete 
with e"^e~ annihilation as a source for f^, i/^, v^^ and i>r neutrinos ("f^"s). The major obstacles 
to obtaining accurate estimates of the emissivity of this process are our poor knowledge of the 
nucleon-nucleon potential, of the degree of suitability of the Born Approximation, and of the 
magnitude of many-body effects (Hannestad and Raffelt 1998; Raffelt and Seckel 1991; Brinkman 
and Turner 1988). Since the nucleons in protoneutron star atmospheres are not degenerate, we 
present here a calculation of the total and differential emissivities of this process in that limit and 
assume a one-pion exchange (OPE) potential model to calculate the nuclear matrix element. To 
acknowledge ignorance, we encourage that a fudge factor of order unity, but perhaps as low as 
0.1, be appended to the rate. The formalism we employ has been heavily influenced by those of 
Brinkman and Turner (1988) and Hannestad and Raffelt (1998), to which the reader is referred 
for details and further explanations. 
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Our focus is on obtaining a useful single-neutrino final-state emission (source) spectrum, as 
well as a final-state pair energy spectrum and the total emission rate. For this, we start with 
Fermi's Golden Rule for the total rate per neutrino species: 

= p')V[n||] p^p^- -^^)- (") 

where 5^(P) is four-momentum conservation delta function, oj is the energy of the final-state 
neutrino pair, {uJu,qu) and {uy^qp) are the energy and momentum of the neutrino and anti-neutrino, 
respectively, and pi is the momentum of nucleon i. Final-state neutrino and anti-neutrino blocking 
have been dropped. 

The necessary ingredients for the integration of eq. (^l]) are the matrix element for the 
interaction and a workable procedure for handling the phase space terms, constrained by the 
conservation laws. We follow Brinkmann and Turner (1988) for both of these elements. In 
particular, we assume for the n + n^n + n + vu process that the matrix element is: 



= -rG if m^) ffA ( ,2 I 2 ) +••• — — = ^ — 2-' (^2) 
4 L fc^ + rn± J 

where the 4 in the denominator accounts for the spin average for identical nucleons, G is the weak 
coupling constant, / (~ 1.0) is the pion-nucleon coupling constant, qa is the axial-vector coupling 
constant, the term in brackets is from the OPE propagator plus exchange and cross terms, k is 
the nucleon momemtum transfer, and m-,^ is the pion mass. In eq. 
terms from the weak part of the total matrix element. To further simplify the calculation, we 
set the "propagator" term equal to a constant C,, a number of order unity, and absorb into C, all 
interaction ambiguities. The constant A in eq. ( ^2|) remains. 

Inserting a j 5{uj — ujy — u}u)doj hy the neutrino phase space terms times bJUJuOJu/^"^ and 
integrating over ujp yields: 



42), we have dropped ■ k 



UJuUjp d^qu (fqp 1 Loljuj-uJyf ^ 

^ ^ /n N^n N-^o ' T^T^ / duJudu , (43) 



u;2 (27r)32w^ (27r)32u;p {2it)^ Jo Jo ^ 
where again uj equals {io^ + Wp). If we integrate over LOy, we can derive the uj spectrum. A further 
integration over uj will result in the total volumetric energy emission rate. If we delay such an 
integration, after the nucleon phase space sector has been reduced to a function of u and if we 
multiply eq. (|4l|) and/or eq. (43) by /a;, an integration over u) from ooy to infinity will leave the 
emission spectrum for the single final-state neutrino. This is of central use in multi-energy group 
transport calculations and with this differential emissivity and Kirchhoff 's Law (§^) we can derive 
an absorptive opacity. 



Whatever our final goal, we need to reduce the nucleon phase space integrals and to do this 
we use the coordinates and approach of Brinkmann and Turner (1988). We define new momenta: 
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p+ = (pi + P2)/2, P- = {pi — P2)/2, P3c = P3 — P+, and p^c = Pa — P+-, where nucleons 1 and 2 are 
in the initial state. Useful direction cosines are 71 = p+ 'P-Zb+lb-l and 7c = p+ ■ P3c/\p+\\P3c\- 
Defining ui = pf 121111 and using energy and momentum conservation, we can show that: 

d?pid?p2 = M^p+d?p- 

UJ = 2T{U- - U3c) 
n3,4 = U++U3c±2{u+U3c)^^^Jc- (44) 



In the non-degenerate limit, the .7-'i^2(l— -^s)!! — -^4) term reduces to e'^'^e~'^^'^+^'^~\ where y is 
the nucleon degeneracy factor. Using eq. (0), we see that the quantity (n+ + n_) is independent of 
both 71 and jc- This is a great simplification and makes the angle integrations trivial. Annihilating 
d^Pi with the momentum delta function in eq. (^), noting that pfdp = QU}!^ — v^J'^dui, pairing 
the remaining energy delta function with and integrating u+ from to 00, we obtain: 

dQnh = 28 X 3 X 5^8.5 ^'•'^''^""^^(^/^)" [ + xuj/Tf'^dx] du . (45) 

The variable x over which we are integrating in eq. (^) is equal to 2u3c- That integral is analytic 
and yields: 

/■oo 

e-^(x2 + xu^/Tfl^dx = rjeiRiirj) , (46) 







where Ki is the standard modified Bessel function of imaginary argument, related to the Hankel 
functions, and rj = uj/2T. Hence, the uj spectrum is given by: 



dQnb -u}/2T i /r,rr\ 

— — oc e ' Lu Ki{uj/2T). 
dio 



(47) 



It can easily be shown that (w) = 4.364T (Raffelt and Seckel 1991). Integrating eq. ( ^5|) over oj 
and using the thermodynamic identity in the non-degenerate limit: 



/ 27r \3/2 



\ 3/2 

n„/2 , (48) 



where is the density of neutrons (in this case), we derive for the total neutron-neutron 
bremsstrahlung emissivity of a single neutrino pair: 

Qnb = 1.04 X l030^(X„pi4)2( J_)5.5 ergs s'^ , (49) 

MeV 

where pi4 is the mass density in units of lO^'' gm cm"'^ and Xn is the neutron mass fraction. 
Interestingly, this is within 30% of the result in Suzuki (1993), even though he has substituted, 
without much justification, (1 + uj/2T) for the integral in eq. (|45|). ([1 + (7rr//2)^/2] is a better 
approximation.) The proton-proton and neutron-proton processes can be handled similarly and 
the total bremsstrahlung rate is then obtained by substituting X"^ + X^ + for in eq. 
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(Brinkmann and Turner 1988). At = 0.7, Xp = 0.3, p = 10^^ cm.'^ , and T = 10 MeV, 
and taking the ratio of augmented eq. ( ^9| ) to the total rate for e+e" production of f^z^^ pairs 
(Dicus 1972), we obtain the promising ratio of ~ bQ. Setting the correction factor Q equal to ~ 0.5 
(Hannestad and Raffelt 1998), we find that near and just deeper than the neutrinosphere, the 
bremsstrahlung rate is larger than that for classical pair production. 



If in eq. (|4^) we do not integrate over LVy, but at the end of the calculation we integrate over 
io from ujy to cxd, after some manipulation we obtain the single neutrino emissivity spectrum: 

dQ'nb .r.(Qnh\ 3 r ^~\. , ^. b^2, .^fQnb\ 3 /"^ e'^^^^^ 



2C{^)u;l £ '—K,mn - nlfdv = 2C{^yi f '-^{e - if'di , (50) 



where rfl = u)y/2T, C is the normalization constant equal to '^^^^ii^^^ (— 0.564), and for the 
second expression we have used the integral representation of Ki{rj) and reversed the order of 
integration. In eq. (|50[), Qnb is the emissivity for the pair. 



Eq. (|50| ) is the approximate neutrino emission spectrum due to nucleon-nucleon 
bremsstrahlung. A useful fit to eq. (|50|), good to better than 3% over the full range of important 
values of r/j^, is: 

dQ'^^ ^ 0.234QnbfiOu\^-\ 



■(^) e---/-. (51) 



dwi, T ^ T 

Setting C equal to 0.5, we have incorporated bremsstrahlung into our Feautrier transport algorithm. 
In j |lO| , we show how the emergent spectrum depends upon (. 



8. Basic Neutrino Transport Results 

The formalism and microphysics described in through §0 were used to calculate the 
neutrino radiation fields for two snapshot profiles in temperature, density, electron fraction, and 
velocity. One of these is from the work of Burrows, Hayes, and Pryxell (1995) and represents the 
wind phase that follows explosion (Model W). The second profile (our Model BM) was kindly 
provided to us by Tony Mezzacappa and is from a multi-group, flux-limited diffusion simulation 
by Bruenn (Messer et al. 1998), 106 milliseconds after the bounce of the core of the Weaver and 
Woosley (1995) 15 Mq star. Since Messer et al. (1998) have already published their results for 
this model, in order to facilitate comparison we highlight our results for Model BM. Note that our 
focus is on neutrino atmospheres and not on completely self-consistent profiles and their evolution. 
Hence, differences between the equations of state and microphysics employed in two different 
dynamical calculations, in particular any differences between the fis, will translate at a given 
epoch into differences in composition and thermal profiles. Post-processing one group's snapshots 
with the code of another can lead to differences in the neutrino fields that are larger than the 
differences in their thermal profiles. The and neutrino luminosity profiles and spectra are 
particularly sensitive to differences between the /Is used. To check this, after achieving a steady 
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state we turned on the Ye coupling for about 5 milliseconds. The upshot was that Ye changed very 
little, demonstrating that we were using substantially the same /is as Messer et al. (1998). 

We concentrate on the generic features of the energy, angle, and spatial distributions of the 
various neutrino radiation fields. We use 50 energy groups, concentrating them below 50 MeV, 
so that the emergent spectra arc wcll-rcsolvcd. The models have 120 spatial grid points out to 
a radius of about 300 kilometers and we interpolate in the various original models to resolve 
important features, such as the neutrinospheres and the shock wave (for Model BM). Since we 
are using the tangent-ray method to set up and determine the angular grid, we employ from 119 
to a few angular groups. In the code, we can establish an arbitrary number of "core rays" in the 
interior to increase the angular coverage at small radii, but we found that we did not need to 
exercise this option. 

The temperature (T), density (p), and Ye profiles for the two models are shown in Figure 1. 
Model BM is a pre-explosion protoncutron star in a stalled shock configuration, while Model W 
is a snapshot of a post-explosion neutrino-driven wind that expands off of the protoneutron star 
after explosion. In Model W from Burrows, Hayes, and Pryxell (1995), Ye asymptotes to a value 
determined by the the competition between Ve neutrino absorption, e~ and capture on 

nucleons, and the speed of expansion. This situation is similar to that found in a gas-dynamic 
laser or freeze out in the early universe. The actual asymptotic Ye and acceleration timescale 
will depend, in a self-consistent calculation, on the details of the neutrino-matter coupling and 
radiation fields and will be the subject of a future paper. Also shown on the lower panel of Figure 
1 are the neutron, proton, and alpha particle mass fractions that bear on the physics of wind 
acceleration and the viability of this wind as a site for the r-process (Woosley and Hoffman 1992; 
Qian and Woosley 1996). 

8.1. Optical Depths and Scattering Albedos versus Radius and Energy 

The integrated depth versus radius or interior mass provides a measure of the global context 
of any transport problem. Figure 2 shows the depth versus radius and neutrino energy for Ve 
neutrinos with energies from 5 to 30 MeV in Model BM. This is not the Rosseland mean which, 
due to the much higher average neutrino energies in the deep interior, reaches a value greater 
than 10^ at the center. The position of the shock wave is manifest. Figure 2 demonstrates 
that the position of the neutrinosphere (r ~ 2/3) is a stiff function of neutrino energy. For i/g 
neutrinos and the energies depicted in the figure, the radii of the neutrinospheres range from ^ 50 
kilometers to ~130 kilometers, more than a factor of two. For the Ve and neutrinos, the range 
is similarly broad, though due to the weaker neutrino-matter coupling for these neutrinos the 
radii are correspondingly smaller. These facts emphasize the dubious merit of referring to a single 
neutrinosphere for a given species. Figure 3 depicts the positions of the neutrinospheres versus 
energy and type. In Model BM, while 10 MeV Ve neutrinos decouple at ~80 kilometers, 100 MeV 
neutrinos decouple as far out as the position of the shock. This situation has a bearing on the 
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strength of the high-energy spectral tail. Note that for the z/e and neutrinos the gain region for 
Model BM, between ~110 kilometers and the shock, resides at optical depths below ~ 0.1 near 
the peak of their respective emergent spectra. For shghtly higher neutrino energies, the optical 
depth of this region is correspondingly higher. Hence, energy deposition in this semi-transparent 
region is problematic and requires a full transport code to study adequately. 

It is important to distinguish absorption from scattering. The scattering albedo is the a priori 
probability that an interaction is a scattering {i^u/xu)- It is a function of composition, neutrino 
energy, neutrino type, and final-state blocking. For v^. neutrinos, the excess of neutrons over 
protons in the free-nucleon, high-entropy region interior to the shock results in an albedo near 
0.25, while for the Ve neutrinos it is 0.5 — 0.6. Figure 4 depicts the Model BM scattering albedos 
versus radius as a function of energy for v^, and v^, neutrinos. In the interior, the absorption 
process, + n — > e~ + p, is suppressed by blocking due to final-state electrons. This results in 
an elevated scattering albedo for the lower energy neutrinos. For neutrinos in Model BM, 
scattering predominates and exterior to 20 kilometers the albedo is above 0.95. Such a scattering 
albedo for the t'^ neutrinos makes its transport a thermalization depth problem that can not be 
easily handled with flux limiters. 



8.2. Emergent Spectra, Luminosities, and Heating Rates 

The emergent neutrino spectra and luminosities are functions of progenitor and they evolve. 
Generally, the spectra after bounce harden with time (Mayle, Wilson, and Schramm 1987), but 
after hundreds of milliseconds or as accretion reverses into explosion (or otherwise subsides), the 
spectra start to soften. The residue then cools inexorably over many seconds, like a clinker plucked 
from a smelter (Burrows and Lattimer 1986). Our Models are merely snapshots, but they serve 
as contexts in which to study the influence of various effects and physics. In addition, the results 
can serve as benchmarks against which to compare those from approximate schemes (see §^). The 
luminosity profiles and spectra for Model BM are depicted in Figures 5 and 6, respectively. The 
Vfj, neutrino luminosity includes that due to f^, v^, Vt, and Ur neutrinos. The steeper rise and 
plateau of the luminosity is a consequence of the small scattering albedo and deeper point of 
energy decoupling, even though the r = 2/3 surface is at larger radii. The peaks in the Ve and Ue 
luminosities mark the inner radius of the gain region, which resides where the luminosity slope is 
negative. The rapid variation in luminosity at smaller radii is a consequence of the variation in 
the temperature slope in the original model, itself presumably a consequence of sparse zoning. 

The asymptotic I'e and v'e luminosities are 4.3 x 10^^ ergs s~^ and 3.1 x 10^^ ergs s~^, 
respectively, 13% higher and 9% lower than the corresponding "BOLTZTRAN" numbers from 
Messer et al. (1998). The differences must stem from a combination of differences in our numerical 
algorithms, in our spatial, angular, and energy zoning, and in our cross sections. Our emergent 
spectra for Model BM are given in Figure 6. The hardness hierarchy of Vf, < i>e < t^^i is manifest, 
as is the dominance of neutrinos at high energies. The i^e and i>e spectra can be fit by a 
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Fermi-Dirac distribution with temperatures and 77s of 2.22 MeV and 3.16 for the Vf. neutrinos and 
2.80 MeV and 3.48 for the neutrinos. The best Fermi-Dirac fit to the z/^ neutrino spectrum has 
a negative ry, which might as well be —00. Note that the emergent v^^ spectrum shown in Figure 6 
was calculated with the bremsstrahlung C set equal to 0.5. The dependence of the spectrum on 
C will be explored in § p!o| . 

The corresponding energy-integrated inverse flux factors (/ Jydej J Hyde) for Model BM are 
plotted versus radius in Figure 7. Figure 8 depicts the unintegrated Ue inverse flux factors {Jy/H^) 
at given radii versus neutrino energy. Since neutrino-matter heating terms are proportional 
to Jy, the higher the inverse flux factor the more efficiently a given energy flux (luminosity) 
heats the matter in the semi-transparent gain region. Different transport algorithms result in 
different inverse flux factors, so getting this term right can be important to the viability of the 
neutrino-driven supernova mechanism (Mezzacappa et al. 1998) and to the acceleration and 
entropy of the post-explosion wind (Burrows 1998a, b). In addition, the harder the spectrum, the 
stronger the neutrino-matter coupling, so the V(, and Ve neutrino spectra versus radius around and 
exterior to the neutrinospheres have a direct bearing on the heating rate. 

Figure 9 portrays the and Jy spectra as the fg neutrinos decouple. As this figure shows, 
at large radii Hy and Jy are the same, but at depth Jy is much larger than Hy. The precise 
values of Jy as the neutrinos decouple determine the matter heating rate. The energy-integrated 
heating and cooling rates versus radius for Model BM for all neutrino species individually are 
given in Figure 10. The positions of radiative equilibrium are indicated with a large dot and the 
inner radius of the gain region for each neutrino is denoted by an X. Note that the gain region 
identified on Figure 10 coincides with the gain region determined from the luminosity plot (Figure 
5). Also included on Figure 10 are the heating rates due to v^, — Vf, annihilation and to — and 
i/t — Vt annihilation, done properly with the appropriate angular factors (Janka 1991). Aside from 
being competitive in the irrelevant unshocked regime, heating due to neutrino pair annihilation 
is meager, at best. In addition, due to the fuzziness of the neutrinospheres, the heating rate per 
cm~^ does not follow the 1/r^ law that might have been appropriate for a sharp neutrinosphere. 
The difference between the heating and cooling curves, the "net gain," for Model BM is given by 
a solid line in Figure 11, to be compared with Figure 8 of Messer et al. (1998). We obtain similar 
heating rates throughout most of the supernova atmosphere, but slightly greater rates between 
110 and 130 kilometers. This slight difference could be due to a combination of things, including 
different techniques, different cross sections, slightly different equations of state, or our better 
angular and energy resolution. 

8.3. Consequences of Various Physical Terms 

The neutrino radiation fields depend upon terms that incorporate various physical effects. 
It is conceptually useful to gauge these terms by their influence on the emergent spectra and 
on the heating rate. Examples of effects that may or may not be included in simpler schemes 
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are the final-state electron blocking term for the charged-current absorption process ( |6.1| ), 
stimulated absorption corrections (§||), weak magnetism and recoil (§|6.1|), and the velocity 
advection, aberration, and Doppler shift terms in the transport equation (eqs. |^ and ^). The 
net gain and the Vf. and neutrino spectra for Model BM, with and without the blocking, 
weak magnetism and recoil, or the stimulated absorption terms, are depicted in Figures 11 and 
12. The blocking correction to the emergent luminosity is ~ 15% and that due to stimulated 
absorption is ^ —3.5%. The blocking and stimulated absorption corrections to the emergent 
neutrino luminosity are of opposite sign and approximately equal to —6.0% and 3.5%, respectively. 
Blocking and stimulated absorption shift the emergent fe and spectra in opposite directions 
in a given energy group by as much as ~20% and ~ —8%, respectively, due to blocking and 
~ —5.5% and ~ 6.5%, respectively, due to stimulated absorption. Blocking increases the net gain 
by 10-20%, while stimulated absorption decreases it by less than 5%. Without electron blocking, 
the absorption cross sections are artifically enhanced. Since the degree of degeneracy is different 
in the envelope and around the neutrinospheres, the magnitude of this effect in the two regions is 
different, enhancing the emergent luminosity more than it decreases the coupling in the periphery. 
Stimulated absorption has the opposite effect (§^), but its effect in the core and in the envelope is 
similarly differential. 

In these Model BM calculations, the effects of weak magnetism and recoil on the emergent Vf. 
and Vf, neutrino spectra and luminosities are small (< 2.0%). This is due in part to the fact that 
the presence of scattering mutes the effect of changes in the absorption cross section through the 
thermalization depth effect. Due to the modest scattering albedo (Figure 4), the response of the 
radial dependence of the radiation field to changes in the absorption cross section is not linear 
with the change in the absorption cross section itself. The increase in the net gain that one would 
anticipate due to any increase in the Vf, luminosity is countered by the concommitant decrease in 
the absorption cross section in the gain region. 

The winds that emerge from protoneutron stars after their envelopes supernova are powered 
by neutrino energy deposition in the expanding gas. Just as with the supernova itself, the wind 
mass and enthalpy fluxes, velocities, entropies, and compositions are influenced by details of 
neutrino-matter coupling and neutrino transport. The distribution of the heating determines the 
magnitude, spatial extent, and timescale of acceleration. In turn, the degree of r~processing in 
the ejecta is a function of the expansion timescale, the asymptotic Y^^ and entropy (Woosley and 
Hoffman 1992; Qian and Woosley 1996). Hence, it is important to gauge the relative strengths of 
the various terms that determine the degree and distribution of neutrino-matter heating. 

The effect of the velocity terms on the emergent Vf, and Vf. spectra for Model W is depicted in 
Figure 13. Model W is the post-explosion wind model from Burrows, Hayes, and Fryxell (1995), 
in which the velocities at large radii are ~30,000kms^^. Of course, at small radii they are zero. 
As Figure 13 shows, the velocity effects collectively boost the emergent spectra of Model W in a 
given energy group by ~15%, with a corresponding boost in the v^. and Vf. luminosities by 15% and 
13%, respectively. This is mostly a consequence of the Doppler shift of the radiation field. Due to 
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the smaller velocities in the important accretion regions in Model BM, the velocity corrections for 
that model are much smaller (< 5%). Figure 14 shows the net gain in Model W for our fiducial 
model, as well as without blocking, velocity corrections, or weak magnetism/recoil corrections 
and implies that various terms not easily or often included in flux-limited or energy-integrated 
transport can each make a ~ 10% difference in the parameters of the wind. Note that though the 
weak magnetism/recoil corrections for Model BM are small, those for Model W are modest. This 
result implies that the importance of absorption cross section changes is a function of the specific 
thermal and composition profiles. What distinguishes the wind is the more abrupt transition 
from the diffusive to the streaming regime and the lower value in its decoupling region (Figure 
1). Whereas, in Model BM the slight increase in the core luminosity due to the inclusion of the 
weak magnetism/recoil term is nullified by the decrease in the absorption opacity in the envelope 
when determining the net gain, in Model W the slight increase in the luminosity due to the 
lower absorption cross section that is a consequence of weak magnetism (particularly for the z^g 
neutrinos) is not adequate to counter the resulting greater transparency of the envelope. The 
Yg at the base of the wind's atmosphere is smaller than that near the neutrinospheres in Model 
BM, with the result that the scattering albedo for the neutrinos is larger there. This results 
in a smaller increase in the emergent luminosity that can't compensate for the increase in the 
transparency of the wind's mantle. 

The anisotropy of neutrino-nucleon scattering and the difference between the transport and 
the total cross sections (eq. |l^) can in principle translate into larger inverse flux factors and, 
hence, greater net gain. Backscattering increases Ji, for a given Hu and delays the transition from 
isotropic to forward-peaked radiation fields. However, since absorption plays an important role 
for the fe and neutrinos and their scattering albedos are not very close to one, the backscatter 
effect is muted. The upshot is that anisotropy accounts for only ~ 2% of the net gain and results 
in shifts of less than 1% in the i^e or z^e spectra. 



9. Flux Limiters 

It is common to seek methods for solving eqs. (12) and (^) without employing the full 
machinery of transport. In principle, such methods simplify the mathematics and speed solution, 
but they compromise accuracy. Eqs. (^) and ([l^) are the first two equations in a moment 
hierarchy in which each moment equation involves still higher-order moments; to solve such 
a hierarchy precisely requires the solution of an infinite number of moment equations. The 
simplifying ansatz often introduced is that higher-order moments can be written in terms of 
lower-order moments, thereby closing the system of equations. However, these so-called closure 
relations can vary a great deal. The most common closure is the flux-limiter. 



In flux- limited schemes, eq. (13) is reduced to its diffusion form in which H,y is set equal to 
the product of the gradient of Pi, and a coefficient (the flux limiter). Pi, is set equal to 1/3 Ji, (the 
Eddington closure), and this expression for H^, is inserted into the divergence term in eq. (|T^). 
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Only this zeroth-moment (energy) equation is solved. The third-order moment, A'^i^, important for 
multi-energy group calculations in moving media, is generally ignored. The art of this approach is 
in the choice of the flux limiter, so-called because the coefficient is constructed in such a way that 
the flux, otherwise mathematically diffusive in this scheme, does not exceed the streaming limit, 
Jy. This is necessitated by the dropping of eq. ([T^), with its causality-enforcing time derivative. 
All angular information is lost, though some flux limiters are derived under certain assumptions 
about the angular distribution (c/. Levermore and Pomraning 1981). Examples of flux-limiters 
that have been employed in supernova calculations are Wilson's (Bowers and Wilson 1982) and 
Bruenn's (1985; Messer et al. 1998). Their prescription for the flux is: 

where Ep = 1 + 3/(1 + \R\/2 + for Wilson and Ep = 1.0 for Bruenn, R = -K^^, 

and is the transport mean-free-path. The expression in parentheses is the limiter. Hence, in 
flux-limiter closures, is a function of only Jy and its derivative (or and R) and (= Pu/ Ju) 
is set equal to 1/3. Note that the Knudsen parameter (R) is small in the diffusion limit. The 
subscript is a reminder that these equations apply for each energy group. It is important to note 
what should be obvious: not all flux limiters are the same, nor are their quantitative consequences. 
Hence, there is not a uniform "flux limiter" result. For instance, the MGFLD scheme of Messer 
et al. (1998) is merely one such approach. The results of employing a given flux limiter can deviate 
from the correct result as much as can the results derived using various flux limiters deviate from 
one another. 



9.1. Angular Distributions 

To gauge the character of the proper angular dependence of the neutrino distribution 
functions in the supernova context, we provide in Figure 15 the Model BM Eddington factor (p^) 
versus radius for z^e neutrinos, at particle energies from 5 to 30 MeV. As we would expect from the 
decoupling hierarchy, the Eddington factors start their rise from the isotropic value of 1/3 from 
the deepest layers. However, for all neutrino species, particularly for the and z/^ neutrinos, the 
Eddington factor is a stiff function of energy and only gradually makes the transition from 1 /3 to 
0.75 over a region that can be 50 to 150 kilometers wide. Many flux limiters effect the related 
transition from diffusion to free streaming early (at higher t^s) and within an unphysically narrow 
range in radius and t^. This can be seen in Figure 16 where the Bruenn and Wilson flux limiters 
for z/e neutrinos at 20 MeV in Model BM are compared with the "effective" flux limiter derived 
using full transport. Approximately 20 kilometers interior to the appropriate point, the Bruenn 
and Wilson limiters begin to deviate from the diffusion value of 1.0. 

Polar plots depicting representative angular distributions of the Model BM z/g specific 
intensity {lu) field for an energy of 15 MeV are presented as solid lines in Figure 17. The transition 



-26- 



from isotropy to forward-peaked is clear, as is the gradual nature of that transition. There is no 
corresponding angular function for either the Bruenn or Wilson limiter. 

Complementary to this polar plot are Figures 18 through 20 of the full transport phase-space 
densities (J'u) versus energy at various radii and for all the neutrino species. Depicted are 
the phase-space densities along the 0° and 90° directions. For the i/g neutrinos, the degree of 
degeneracy at depth is clear; one can almost read the Vg chemical potentials off the graph. From 
Figure 19, we see that the occupancy of the Ug neutrino states is generically low, but from Figure 
20 we see that at depth and for low energies the occupancy of the neutrino states can approach 
0.5. Blocking due to final-state occupancy is generally unimportant in the pair source terms, 
since the peak energies of the pair source functions are always significantly above the energies at 
which J^^ is high. 

Among other things, Figures 18 through 20 convey a sense of the angular dependence of J-'u, 
ignored in the standard flux-limitcr schemes. At depth, since the radiation fields arc isotropic, the 
0° and 90° curves are the same. However, with increasing radius and at lower energies, deviations 
from isotropy manifest themselves; transverse beams are less occupied than forward beams. As 
expected, at low optical depths this differential effect is quite pronounced. Flux-limited transport 
schemes are not capable of addressing or illuminating this phenomenology. 

9.2. Neutrino Heating and Emergent Spectra Using Flux Limiters 

Given the nuances in the angular distributions of the neutrino radiation fields portrayed in 
Figures 15-20, it is no wonder that flux limiters only inadequately represent the radiation energy 
density profiles, emergent spectra, and net gain in the semi-transparent decoupling region. This is 
made manifest by comparing the emergent spectra and net gain derived using such flux limiters 
with those same quantities obtained using full transport. Figure 21 compares the emergent i/g 
spectrum for the BM model using the full Feautrier /tangent-ray formalism, Bruenn's limiter, and 
Wilson's limiter. Though comparisons for many snapshot profiles would be useful, we can still 
conclude for such an early protoneutron star epoch that Bruenn's limiter, as simple as it is, results 
in a spectrum that deviates from the more precise spectrum by '^5-10%, while Wilson's limiter, 
despite its modestly greater complexity, can be off by as much as ~20-30%. Figure 22, in which 
curves of the net gain (heating) versus radius are compared, tells a similar story: Bruenn's limiter 
yields net gains that are generally off, but by no more than ~20%, while Wilson's can be off by as 
much as ~50%. Figures 21 and 22 serve to illustrate both that all fiux limiters are not the same 
and, in particular, that they can underestimate the net gain in the outer gain region by many tens 
of percent. 

In sum, flux-limiter schemes can miscalculate net heating rates, radiation energy densities 
(quantities that factor into the net gain) , emergent spectra, and the inverse flux factors by from 5% 
to 50% and can artificially accelerate the transition from isotropy to free-streaming in the r < 1 
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region. This is particularly true for neutrinos, with their extended neutrinospheres. Furthermore, 
the thermalization depth effect is difficult to handle with flux limiters when the scattering albedo 
is large. The albedo for neutrinos is above 0.90 throughout most of the object. As a result, 
only full transport can properly handle the enhancement in the effective absorption path due to 
the frustrated escape caused by scattering. 



10. Determinants of the Emergent Neutrino Spectrum 

The Ufj, and Vr neutrinos and their antiparticles carry away from the protoneutron star 
more than 50% of its total binding energy. Since they do not participate in charged-current 
interactions, they energetically decouple at smaller radii and, hence, at larger temperatures, than 
the other neutrino species. This results in a harder spectrum (Figure 6) and the hardness hierarchy 



alluded to in §8^. The fact that there are four species is primarily responsible for their major 
cooling role. Neutrino-matter energy coupling is affected by the inverse production processes 
of pair annihilation and nucleon-nucleon bremsstrahlung (§0), as well as by neutrino-nucleon 
and neutrino-electron scattering. The proper treatment of energy redistribution by scattering 
is deferred to a later publication. However, it is clear that scattering generically softens the i/^ 
spectra. 

Ignoring the potential effects of neutrino oscillations, the emergent spectra have a direct 
bearing on the process of neutrino nucleosynthesis (Woosley et al. 1990) and on the suitability 
of various underground detectors that rely on neutral-current spallation processes with high 
energy thresholds. In both cases, the relevant neutral-current interaction cross sections are stiffly 
increasing functions of neutrino energy, with thresholds above ~ 15 MeV (Haxton 1990). Hence, 
they are most sensitive to the i/^ component and its precise spectrum. In the past, people had 
thought that the spectra were hard, with effective Fermi-Dirac temperatures of ~ 8 — 9 MeV 
and average energies of ~ 25 — 30 MeV. However, the energy spectrum on Figure 6 can be very 
approximately fit with a temperature of 7 MeV. 

Bremsstrahlung has a major effect on the radiation field. The factor ^ in incorporates 
a correction for our approximations to the propagator terms and to the nuclear matrix element. 
In Figure 6, C was set equal to 0.5. Using Hannestadt and Raffelt (1998) and our own estimates 
of the correct propagator terms, we derive that above 10^^ g cm~^ C, is above 0.7 and that at and 
around 10^^ g cm^^ ^ is near 0.2. This translates into an "average" C of ~ 0.5 for protoneutron 
stars. Figure 23 depicts the consequences of varying from 0.0 to 1.0 in steps of 0.2 for the 
emergent energy spectrum. Due to the presence of an absorption term for every emission term 



(Kirchhoff's Law; eq. 36), the strength of the spectrum is not strictly linear in As Figure 23 
demonstrates, nucleon-nucleon bremsstrahlung is softer than e~^e~ annihilation (the other major 
u^Ufi source) and can dominate at low energies. Though the emergent spectra are softer, due to the 
extra source the spectra are also brighter at every energy. Hence, the inclusion of nucleon-nucleon 
bremsstrahlung increases the flux, while decreasing the average and peak neutrino energies. This 
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is important. At 10 MeV, the spectrum can be more than a factor of two stronger with 
bremsstrahlung than without. For energies above ~ 35 MeV, e^e~ annihilation stih dominates 
the emergent spectrum. In Figure 23, the lowest curve corresponds to a pure e+e" annihilation 
source. Note that it is demonstrably harder than when C, is large and that it alone is "pinched." 
Though it still remains to be determined whether nucleon-nucleon bremsstrahlung in supernovae 
is in fact dominant for spectrum formation, Figure 23 suggests that it is, particularly at lower 
neutrino energies. Since energy transfer due to neutrino— matter scattering and gravitational 
redshifts will only further soften the emergent spectra, we conclude that spectra are indeed 
softer than traditionally quoted. 

Also shown on Figure 23 is an emergent spectrum with the scattering cross sections very 
artificially cut by one half. This curve demonstrates the severe dependence of the spectra on the 
basic numbers associated with the neutrino-matter interaction. It suggests that if we did not have 
a fairly good handle on the basic interactions of neutrinos with nucleons our predictions would be 
quite different, and perhaps would be quite wrong. 

11. Summary 

We have constructed and described an implicit, multi-group, multi-angle, multi-species 
neutrino transfer code to be used in the context of core-collapse supernovae and protoneutron 
stars. The basic algorithm embodies the Feautrier and tangent-ray approaches to spherical 
atmospheres and is conceptually equivalent to various Boltzmann solvers. It is capable of 
resolving angular distributions and of calculating angular moments to great precision and employs 
accelerated A iteration to achieve rapid convergence. Focusing on "neutrino atmospheres," we 
presented the energy spectra, neutrino heating rates, Eddington factors, angular distributions, 
and phase space densities for typical post-bounce structures. The influence on these quantities, 
in particular on the net gain, of various corrections to the charged-current cross section and 
terms in the transport equation were examined and the character of the neutrino radiation fields 
and spectra was scrutinized. One goal has been to provide a detailed snapshot of the neutrino 
radial, angular, energy, and species distributions in a typical post-bounce environment, including 
in the protoneutron star context, and to explore the factors that determine the heating rates in 
the semi— transparent gain region, so central to the viability of the neutrino-driven mechanism 
of supernova explosions. To this end, we focused on the decoupling transition of the emergent 
neutrinos. Moreover, we compared the emergent spectra and neutrino heating rates obtained 
using representative flux limiters with those obtained using our Feautrier transport algorithm 
to gauge the accuracy of those oft-used approximate schemes. Finally, we derived the rate of 
nucleon-nucleon bremsstrahlung and its neutrino source spectrum and showed for the first time 
that it probably dominates z/^ neutrino production and spectrum formation. 

The tool that we have developed is meant to explore supernova explosions, protoneutron star 
cooling, the neutrino signature of core-collapse, neutrino shock break-out, and post-explosion 
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winds, among other things. It is also easily converted into a photon transport code for the 
study of classical supernova light curves. However, we have yet to generalize the scheme for use 
in multi-dimensional supernova simulations or in the general relativisitic context, nor have we 
parallelized it for use on shared-memory machines. Hence, much technical work remains. 

Supernova theory has been evolving for thirty years and in that time our understanding 
of the neutrino and its interactions has changed substantially. There are now indications from 
atmospheric and solar neutrino experiments that lepton number is not strictly conserved and that 
neutrinos may mix. Heating in the protoneutron star mantle is a subtle sum of competing effects. 
We have investigated in this paper but a few of these. This effort to fully characterize the neutrino 
radiation fields is part of a larger effort, as yet unfinished, to understand the mechanism of 
supernova explosions and their systematics. However, when this puzzle box is eventually opened, 
precise neutrino transport will certainly be one of the keys. 
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Fig. 1. The temperature (T), density (p), electron fraction {¥(,), proton fraction {Yp), neutron 
fraction (Yn), and alpha fraction (Ya) for Models BM (top panel) and W (bottom panel). In Model 
BM, the shock is located at 170 km, but in Model W there is no shock on the grid. 

Fig. 2. — The Ue neutrino optical depth (ti,^) versus radius (in kilometers) for Model BM, at various 
particle energies. As the energy of the neutrino increases the degree of transparency decreases. The 
dip in the optical depth at 170 km is where the shock (and, hence, a density jump) is located. The 
solid horizontal line shows where T,y^ = 2/3. 

Fig. 3. — The neutrinosphcrc radii versus neutrino energy for i/g, u^, and "f^" neutrinos. For 
a given neutrino energy, the 'V^" neutrinos decouple first, resulting in a r = 2/3 radius that is 
smaller than that for either or v'g neutrinos. The hierarchy in decoupling radii of f e > z^e > 
is manifest. 

Fig. 4. — The Ve and Ue scattering albedos versus radius for neutrino energies of 10, 20, and 40 
MeV. The shock is at ~170 km. The increases in the albedo at small radii can be traced to e~ 
blocking of I'e absorption on neutrons, predominantly. 

Fig. 5. — Model BM Vg, Ve, and t'^ luminosities versus radius (in km). The "f^" luminosity is the 
sum of the f^, i^^, i^t, and v-r neutrino luminosities. The modest peaks mark the inner radius of 
the gain region, in which, due to net absorption, the luminosity slope is gently negative. 

Fig. 6. — The Model BM emergent neutrino luminosity spectra for the three neutrino types. The 
symbols indicate the positions of the energy groups (ye - filled squares; Ve - filled triangles; - 
open circles). 

Fig. 7. — The energy-integrated inverse flux factor {{1/Fu)) as a function of radius for i/^ and 
neutrinos. The sharp increase in (1/-Fvef) at 110 km occurs just inside the gain radius, where the 
neutrinos are starting to decouple from the matter. At large radii (off the plot), the inverse flux 
factors approach the unity expected for the free-streaming regime. 

Fig. 8. — The Model BM ratios of the neutrino energy density to the Vg energy flux for radii of 
100, 125, 150, 180, 200, 220, and 300 km. 

Fig. 9. — The Ug neutrino energy flux (Fi,, solid) and energy density (cEu, dashed) spectra at 
various radii. A, B, C, D, and E denote radii at 50, 130, 200, 250, and 300 km. At depth, the 
spectra are very different, but they converge at large distances from the neutrinospheres. At these 
large distances, the unintegrated flux factor, F^/cE^, is unity. 

Fig. 10. — Model BM heating and cooling rates (in ergs g~^ s~^) versus radius (in km). The heating 
and cooling rates for the three neutrino species are shown, along with the u — i> annihilation energy 
deposition rates. The solid points indicate where radiative equilibrium is achieved for each neutrino 
species. The X's indicate the positions of the gain radii for the respective neutrino types. The top 
two solid lines are the heating and cooling curves for the Ug neutrinos. The dashed lines are the 
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heating and cooling curves for the neutrinos. The bottom two sohd hues are the heating and 
coohng curves for neutrinos. The bold dashed curves are the heating rates for the VeVf, — > e^e~ 
process (top) and both the v^v^ — > e+e~ and v^-Vt — > e+e" processes (bottom). 

Fig. 11. — The net heating rate (net gain) for various BM models versus radius. The fiducial 
model (solid) is compared to models with no stimulated absorption (short-dashed line), with no 
e~ blocking (long-dashed line), or with no weak magnetism/recoil (dot-dashed line). The absence 
of either stimulated absorption or weak magnetism/recoil would result in an increase in neutrino 
absorption and, thus, a greater heating rate. The absence of blocking would result in a decrease 
in the net gain. 

Fig. 12. — The emergent luminosity spectrum for both the and v^, neutrinos for our fiducial 
Model BM (solid line), compared with models without stimulated absorption (small-dashed line), 
e~ blocking (long-dashed line), or weak magnetism/recoil (dot~dashed line). Also included is a 
model for which the total scattering cross section is substituted for the transport cross section 
(short dash- long dash). 

Fig. 13. — The emergent v^. luminosity spectrum for the wind Model W with a velocity field (solid 
line) and without a velocity field (dashed line), versus energy in MeV. The velocity terms boost 
the spectrum that emerges from a protoneutron star with a wind by ~15%. 

Fig. 14. — The net heating rate (net gain) versus radius (in km) for wind Model W. The fiducial 
model (solid) is compared with models for which either weak magnetism/recoil (short-dashed line) 
or e~ blocking (long-dashed line) is ignored. Also included is a model for which the velocities were 
set equal to zero (dot-dashed). 

Fig. 15. — The Eddington Factor for v^. neutrinos versus radius (in km) at various neutrino 
energies. At depth, in the diffusive region the Eddington factors converge to 1/3. At large radii, 
the Eddington factors approach unity. The low-energy neutrinos are the first to decouple and their 
Eddington factors approach unity faster than those of the higher-energy neutrinos. 

Fig. 16. — Electron neutrino flux limiter profiles versus radius (in kilometers) for model BM using 
Bruenn's flux limiter (dot-dahsed), Wilson's flux limiter (solid), and an artifical flux limiter derived 



from the full formalism (dashed). The fe neutrino energy is 20 MeV. See text and eq. (52) for 
details. 

Fig. 17. — Polar plots of the specific intensity {ly) of Vf. neutrinos with an energy of 15 MeV. 
Shown are angular distributions at radii of 80, 120, 170, and 300 km. In the interior, the radiation 
fields are isotropic and strong. At large radii, the distribution becomes more forward-peaked and 
geometric dilution decreases ly. The numbers on the left axis provide the scale, with the negative 
numbers emphasizing the fact that the rays are pointing backward in this hemisphere. 

Fig. 18. — The phase-space density [Tu) for the v^, neutrinos versus neutrino energy in MeV, for 
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various radii from 20 km to 170 km. The solid lines are for the forward direction and the dashed 

lines are for the transverse direction (at ~ 90° to the radial direction). At small radii, the Vf. 
neutrinos arc degenerate, but at larger radii, and generally at larger energies, they quickly become 
non-degenerate. Note that at small energies, "larger" radii, and large angles, the degeneracy of the 
Ve neutrinos diminishes. 

Fig. 19. — The phase-space density [Ty) for the i>e neutrinos versus neutrino energy in MeV, at 
radii of 120, 150, and 200 km. The solid lines are for the forward direction and the dashed lines 
are for the transverse direction (at ~ 90°). 

Fig. 20. — The phase-space density {Ty) for the neutrinos versus neutrino energy in MeV, at 
radii of 55, 70, and 80 km. The solid lines are for the forward direction and the dashed lines are 
for the transverse direction (at ~ 90°). At depth, and at lower energies, neutrino degeneracy 
(.Tv) approaches 0.5, as expected for the situation with no net lepton number and, hence, zero 
chemical potential. 

Fig. 21. — The emergent v^, neutrino luminosity spectrum for model BM using the full 
Feautrier /tangent-ray formalism (solid), Bruenn's limiter (short dashed), and Wilson's limiter (long 
dashed) (cf. Figure 6). 

Fig. 22. — A comparison of the net gain (in erg g^^ s^^) versus radius (in kilometers) for model 
profile BM, calculated using the full transport formalism of this paper (solid), Bruenn's flux limiter 
(dot-dashed), and Wilson's flux limiter (dashed) (cf. Figure 11). 

Fig. 23. — The emergent luminosity spectra for Model BM for bremsstrahlung factors, of 0.0, 
0.2, 0.4, 0.6, 0.8, and 1.0. Also included is the spectrum with the f^-nucleon scattering cross 
section artiflcially decreased by 50%. 
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